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Abstract. The solution of the mean spherical approximation (MSA) integral 
equation for isotropic multicomponent dipolar hard sphere fluids without external 
fields is used to construct a density functional theory (DFT), which includes external 
fields, in order to obtain an analytical expression for the external field dependence 
of the magnetization of ferrofluidic mixtures. This DFT is based on a second-order 
Taylor series expansion of the free energy density functional of the anisotropic system 
around the corresponding isotropic MSA reference system. The ensuing results for the 
magnetic properties are in quantitative agreement with our canonical ensemble Monte 
Carlo simulation data presented here. 



1. Introduction 

Ferrofluids are colloidal suspensions of single domain ferromagnetic grains dispersed 
in a solvent. The stabilization of such suspensions is usually obtained by coating the 
magnetic particles with polymer or surfactant layers or by using electric double layer 
formation. Since each particle of a ferrofiuid possesses a permanent magnetic dipole 
moment, upon integrating out the degress of freedom of the solvent, which gives rise to 
effective pair potentials, dispersions of ferrocoUoids can be considered as paradigmatic 
realizations of dipolar liquids fTj. The effective interactions of such magnetic particles are 
often modeled by dipolar hard-sphere (DHS) |2l|3], dipolar Yukawa jl], or Stockmayer 
[5] interaction potentials. The most frequently applied methods to describe ferrofluids 
encompass mean field theories [HI [7], thermodynamical perturbation theory |2], integral 
equation theories [HI El [10], various DFTs [IH [121 [131 [H HSl g] , as well as Monte Carlo 
[T6l [T7] and molecular dynamics [TH [191 120] simulations. 

Within the framework of DFT and the mean spherical approximation (MSA), 
previously we have proposed an analytical equation [1] for the magnetic field dependence 
of the magnetization of one-component ferrofluids, which turned out to be reliable as 
compared with corresponding Monte Carlo (MC) simulation data. For this kind of 
system the effect of an external magnetic field has been taken into account by a DFT 
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method, which approximates the free energy functional of the anisotropic system with 
an external field by a second-order Taylor series expansion around the corresponding 
isotropic reference system without an external field. The expansion coefficients are the 
direct correlation functions which for the studied isotropic dipolar hard-sphere (DHS) 
and dipolar Yukawa reference systems can be obtained analytically from Refs. |211 122] . 

However, in practice the magnetic colloidal suspensions are often multicomponent. 
In order to describe the magnetization of ferrofluidic mixtures we extend our one- 
component theory to multicomponent systems. This extension is based on the 
multicomponent MSA solution obtained by Adelman and Deutch [23]. They showed 
that the properties of equally sized hard spheres with different dipole moments can be 
expressed in terms of those of an effective single component system. Because MSA is a 
linear response theory Adelman and Deutch could predict only the initial slope (or zero- 
field susceptibility) of the magnetization curve. Using their MSA solutions as those of a 
reference system, in the following we present DFT calculations of the full magnetization 
curves of equally sized, dipolar hard sphere mixtures. These results are compared with 
MC simulation data. 



2. Microscopic model and MSA solution 

We consider dipolar hard-sphere (DHS) fluid mixtures which consist of C components. 
The constitutive particles have the same diameter a but the strength rria of the embedded 
point dipole can be different for the components a = 1, C In the following the indices 
a and b refer to the components while the indices i and j refer to individual particles. 
The system is characterized by the following pair potential: 

«;g^^(ri2, ui, U2) = <''(ri2) + <''(ri2, Wi, 002), (1) 

where wj^^ and wj^^ are the hard-sphere and the dipole-dipole interaction pair potential, 
respectively. The hard-sphere pair potential given by 

HSr \ J 00 , ri2 < (T , . 

^'"^ = I , r,2 > a. 
The dipole-dipole pair potential is 

W^^ (ri2,Wl,W2) = 3— L'(Wi2,CUl,CU2), (3) 

^12 

with the rotationally invariant function 

D(a;i2,a;i,a;2) = 3[mi(a;i) ■ r 12] [1112(^^2) ■ ri2] - [mi(c^i) ■ m2(w2)], (4) 

where particle 1 (2) of type i (j) is located at ri (r2) and carries a dipole moment of 
strength nii {rrij) with an orientation given by the unit vector rni{ui) (m2(a;2)) with 
polar angles Ui = {9i, 0i) (a;2 = {O2, ^2)); t'u = ti — r2 is the difference vector between 
the center of particle 1 and the center of particle 2 with ri2 = |ri2|. 

Within the framework of MSA Adelman and Deutch [23] presented an analytical 
solution for the aforementioned C-component isotropic dipolar fluid mixture in the 
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absence of external fields. The importance of their contribution is that it provides 
simple analytic expressions for correlation functions, the dielectric constant (in our case 
the zero- field magnetic susceptibility), and thermodynamic functions. In our envisaged 
DFT calculations for dipolar mixtures with external fields we consider the isotropic 
DHS fluid mixture without external field as a reference system which is described by 
the following MSA second-order direct correlation function: 

c^^\ri2,p,n)D{uji2,uJi,UJ2) + c^l\ri2,p,n)A{uji,uj2) , (5) 
where p = N/V = X]^=i Pa is the total number density in the volume V of the system, 

(2) 

Pa = Na/V is the number density of species a, c\jg is the one-component hard sphere 

I2\ I2\ 

direct correlation function, while c)-, and c\ are correlation functions determined by 
Wertheim [21] for the one-component dipolar MSA fluid at the same temperature but 
evaluated for an effective dipole moment m = m and at an effective packing fraction 
Tj = rj. Accordingly, in Eq. (|5]) we have introduced 



fn=\l ^-=X'\ r^='-a^^^^^y^ (6) 



E a=l ml ^_ ^ ^^ Y.a=im lpa 

C ' 6 rh? 

and the rotationally invariant function 

A(wi,W2) = mi(wi) ■ m2(w2)- (7) 

In order to explain dependence on n it is convenient to introduce a new parameter 
^ = ^(Xl) = which is given by the implicit equation 

4vrx, =g(20-g(-0, (8) 

which has the same form as the corresponding equation for the one-component system 
(see Refs. [HET]) where 

1 ^ 

a=l 

is the averaged Langevin susceptibility and /3 = l/lksT) is the inverse temperature with 
the Boltzmann constant ks- The function q{x) is the reduced inverse compressibility 
function of hard spheres within the Percus-Yevic approximation: 

:i-x) 



<li^) = '-TT^- (10) 



For the zero-field (initial) magnetic susceptibility x of the mixture the theory by 
Adelman and Deutch [23] yields 

X = — ^ — • (11) 

?(-e(xJ) 

Concerning Eq. ^ Wertheim [2T] and Adelman and Deutch [23] showed that 

CaVi2,P,^) = 2n[cgjj(ri2,2np) - c§s{ri2, -np)], (12) 
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CD\ri2, P,n) - 3r^^ dss^c§ {s, p,n) (13) 

^0 



with 

c§{ri2,p,n) = n[2c^^'g{ri2,2np) + c^!;(ri2, -np)], (14) 

where c Hsif 12 ■, p) is the one- component hard sphere Percus-Yevick correlation function 
at density p. The dimensionless quantity n = ^/rj determined by solving Eq. ([8]). We 
find that n vanishes in the nonpolar limit — )■ for all a. This can be inferred from 
the results of Rushbrooke et al. [24j (obtained originally for one-component dipolar 
fluids) according to which the solution of Eq. ^ can be expressed as a power series in 
terms of Xl'- 

^=lx.-'-^xl+Oixl). (15) 
From Eqs. ( ITSi) and ([6]) it follows that 
lim n = 

{mi,...,mc}->-0 

lim i= lim (1 + 0(£J) I Vm^ I =0. (16) 



a=l 



/ON 

Therefore in this limit the functions c}j and vanish (see Eqs. ( |T2|) . ( [T3l) . and ( 1T4|) ). 
Due to < C^^^ < C, as expected in the nonpolar limit the rhs of Eq. 

(2) 

reduces to the direct correlation function c)jg of a one-component HS fluid with total 
density p = Y.a=i P^- 

For a one-component system (C = 1) the prefactor niamij/w? equals 1 and 
rj = 7] = ^pc^ so that Eqs. ([H]) and render ^ for the one-component system which 
indeed yields the one-component direct correlation function. 

Considering the case of a binary mixture of hard spheres and of DHS reveals the 
approximate character of Eq. ([5]). In this case, for the dipolar hard sphere - hard 
sphere cross correlations (i.e., nia 7^ 0, = 0), according to Eq. ([5]) the corresponding 
correlation function reduces to the a one-component hard sphere direct correlation 
function, which certainly is a rough approximation. We note that c^^\ri2 — ?■ 00) ~ 
is long-ranged (see Eq. (fT3l) ) while c^^ (ri2 > ex) = is short-ranged (see Eq. ( fT2l) ). 

In the following we consider DHS mixtures in a homogeneous external magnetic 
field H, the direction of which is taken to coincide with the direction of the z axis. For 
a single dipole the magnetic field gives rise to the following additional contribution to 
the interaction potential: 

u^' = -miH = -TTiiH cos Oi, (17) 

where the angle 6i measures the orientation of the i-th dipole relative to the field 
direction. 
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3. Magnetization in an external field 

In the following we extend our previous theory |i4j to C-component and polydisperse 
dipolar mixtures in which the particles have the same hard sphere diameter but different 
strengths of the dipole moments. 



3.1. Multicomponent systems 

Our analysis is based on the following grand canonical variational functional Q, which 
is an extension to C components of the one-component functional used in Ref . [1] : 

c 

Q = FDHs[Pl,---,PC,{(^li^),---,(^c{^)},T]-^Pa / duaa{u){fla- K'^^U)), (18) 

a=l 

where F^hs is the Helmholtz free energy functional of an anisotropic, dipolar, equally 
sized hard sphere fluid mixture and where Ha and aa{u}) are the chemical potential and 
the orientational distribution function of the species a, respectively. Since the external 
field is spatially constant, pi,...,pc are constant, too. Thus Fdhs is a function of 
pi, pc,mi, ...,mc and a functional of ai{co), ...,ac{oj). The Helmholtz free energy 
functional consists of the ideal gas and the excess contribution: 

Fdhs = i^^lpi, ...,Pc, ttcM}, T] 

+FdhsW PC, n (19) 

For the C-component mixture the ideal gas contribution has the form 

c 



F'^ = kBTVj2pa HpaK)-l+ f dLoaaitu)\n{A7iaaitu)) 

a=l L ^ 



(20) 



where is the de Broglie wavelength of species a. If the system is anisotropic the DHS 
free energy F^^'^^ is approximated by a second-order functional Taylor series, expanded 
around a homogeneous isotropic reference system with bulk densities pi, ...,pc and an 
isotropic free energy F^^'^: 

PPDiif [Pi , • • • , PC, {«! (w) , • • • , ac(^^) } , = ^F^j^h's {pi,-;Pc,T) 
c 



~^P'' j d^riduAaa{u})c'^\pi, pc,T) 

a=l ^ 



^ PaPfe / d^ridui / d^r2du2Aaa{u}i)Aab{u2)c''^il{ri2,u}i,uj2,Pi,---,Pc,T), (21) 

a,b=l ^ ^ 

where Aaa(w) = aa^uj) — l/(47r) is the difference between the anisotropic {H ^ 0) and 
the isotropic (if = 0) orientational distribution function of the component a; ch!^ and 
^ah (s^^ Eq.([5])) are the first- and second-order direct correlation functions, respectively, 
of the components of the isotropic DHS mixtures. Since in the isotropic system all 
Sa^ are independent of the dipole orientation and because \ duiAaa{ui) = 0, only 

(2) 

the second-order direct correlation functions c^^ provide a nonzero contribution to the 
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(2) 

above free energy functional. Since c)j^ depends only on the difference vector ri — r2, 
Eq.( l2Ti) reduces to 

(^^DHs' [pu PC, ■■■,ac{u})},T] = PF^^^'s [pi, pc,T) 

-^pV ^ XaXfo J dui J duj2Aaa{uJi)AabiuJ2) j d^ruc'-^fl {ri2, uji,uj2, Pi, PC, T), 

a, 6=1 

(22) 

where Xa = Na/N = Pa/p is the mole fraction of the component a. The expression for 
the excess free energy PF^^'^ of the isotropic DHS system was also given by Adelman 
and Deutch [23J. In the case of a cylindrical sample (elongated around the magnetic 
field direction ) and homogeneous magnetization all aa{u}) depend only on the polar 
angle 6, and thus they can be expanded in terms of Legendre polynomials: 

1 1 °° 

aa{u;) = —aa{cos9) = —y^ aaiPiicosO), a = 1, 2, C . (23) 

1=0 

Due to aao = 1/2 one has 

^ oo 

Aaa{co) = —y^aaiPi{cose). (24) 

1=1 

The second-order MSA direct correlation functions of the DHS fluid mixture (see Eq. 
([5])) are used to obtain the excess free energy functional. In order to avoid depolarization 
effects due to domain formation, we consider sample shapes of thin cylinders, i.e., needle- 
shaped volumes V. Due to the properties of D, A, and Pi only the terms Oai with / < 1 
contribute to the excess free energy. Elementary calculation leads to 

pexc,ai c^^2 ^ C' 

= foHS - ~ '?(~^ )) 5Z XaXbrnarribaaiabi, (25) 

-^■f- a,b=l 

where f^^g = Fdhs/^- Minimization of the grand canonical functional with respect 
to the orientational distribution functions (note that f^^s does not depend on them) 
yields 

aaico) = Z-^exp {(3ma{H + ^ (1 - g(-0) > >f,m,«w 1 Pi(cos0) 1 , (26) 



(3ma + ^ (1 - q{-^))J2xbmbab?j Pi(cos0) j 



with normalization constants Za which are fixed by the requirements j duaa{uj) = 1. 
With this normalization the expansion coefficients aai are given by 



3, 



R I . 2(l-g(-0) ^ 
finia I H H — 2^ pb^ibUbi 



1,2,...,C, (27) 



where L{x) = coth(x) — 1/x is the Langevin function. Each particle of the magnetic 
fluid carries a dipole moment which will be aligned preferentially in the direction of the 
external field. This gives rise to a magnetization 

c ^ ^ c 

imaaai. (28) 



M = ^ PaTTla j duaa{uj) COS 6* = ^ ^ p^/ 
a=l a=l 
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(29) 



Equations ( 125]) and ( E7|) lead to an implicit equation for the dependence of the 
magnetization on the external field: 

c 

M = p^TUaXaL 
a=l 

We note that in Eq. ( I29l) in the limit of weak fields the series expansion of the Langevin 
function, L{x — )■ 0) = x/3, reduces to Eq. ffTTl) for the zero-field magnetic susceptibility. 

3.2. Polydisperse systems 

For ferromagnetic grains the dipole moment of a particle is given by 

m(x) = ^MsV^ (30) 
o 

where Mg is the bulk saturation magnetization of the core material and V is the 
diameter of the particle. Accordingly, our model of equally sized particles with different 
dipole moments applies to systems composed of materials with distinct saturation 
magnetizations. Another possibility consists of considering particles with a magnetic 
core and a nonmagnetic shell, which allows one to vary m via changing the core size 
with Ms fixed and by keeping the overall diameter of the particles fixed via adjusting the 
thickness of the shell. For a small number C of components this can be experimentally 
realizable. 

Equation (fTT]) has been extended even to the description of polydisperse ferrofluids 
[01 [19]. However, it is unlikely that this extension relates to a realistic experimental 
system because it supposes again that the diameters of all particles are the same. The 
two possible realizations mentioned above will be very difficult to implement for a large 
number C of components, mimicking polydispersity. If one nonetheless wants to study 
such kind of a system the expression for its zero-field susceptibility is a natural extension 
ofEq. 0: 

1 f°° 

X, = ^Pp dVp{V)m\V), (31) 

where p{T>) is the probability distribution function for the magnetic core diameter. The 
corresponding zero-field susceptibility of the polydisperse system is 

X = -^^, (32) 
where ^ is the implicit solution of the equation 

4vrx, =g(20-g(-0- (33) 

We note that similarly the above equation for the magnetization (see Eq. ( l29i) ) can 
also be extended to polydisperse fluids leading to magnetization curves M{H) defined 
implicitly by 



M = p dVp{V)m{V)L f3m(V) (^H + ^^—^^^-^Ad^ 



(34) 
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Figure 1. Zero-field susceptibility x (Eq. (ITlT) ) as function of the averaged Langevin 
susceptibility Xl ^ 1/^ (Eq. In terms of these quantities DFT (MSA) predicts 

the master curve given by the full line. The MC data correspond to the following 
choices of the system parameters. The dipole moments are m\ = 0.5 and rrij = 1 in all 
cases, the reduced densities are p* — 0.6 (diamonds), p* = 0.65 (circles), and p* — 0.7 
(triangles) with the concentrations xi xi — 0.25, xi = 0.5, xi — 0.75, and xi — 1. 
The error bars of the MC data are given by the symbol sizes. 



In the limit of weak fields Eq. (!34ll reduces to the expression in Eq. ( l32ll for the zero- 
field susceptibility. In the following we shall do not assess via MC simulations the range 
of validity of Eq. for polydisperse magnetic fluids, leaving this for future studies. 

4. Monte Carlo simulations 

In order to assess the predictions of the DFT presented in Sec. [3] we have carried out 
MC simulations for DHS fluid mixtures using canonical (NVT) ensembles. Boltzmann 
sampling and periodic boundary conditions with the minimum-image convention [25] 
have been applied. A spherical cutoff of the dipole-dipole interaction potential at half 
of the linear extension of the simulation cell has been applied and the reaction field 
long-ranged correction [25j with a conducting boundary condition has been adopted. 
For obtaining the magnetization data, after 40.000 equilibration cycles 0.8-1.0 million 
production cycles have been used involving 1024 particles. In the simulations with an 
applied field the equilibrium magnetization is obtained from the equation 



where the brackets denote the ensemble average. In simulations without external 
field the zero-field magnetic susceptibility has been obtained from the corresponding 




(35) 
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Figure 2. Magnetization curves of a binary DHS fluid mixture for five values of the 
concentration xi of the species with = 0.5; = 1. The overall number density is 
p* = 0.6. The curves correspond to the predictions of the MSA based DFT (Eq. ([^ '1. 
The symbols are MC simulation data. Their error bars are given by the symbol sizes. 



fluctuation formula 

where Ai = YliLi is the instantaneous magnetic dipole moment of the system. 
Statistical errors have been determined from the standard deviations of subaverages 
encompassing 100.000 MC cycles. 



5. Numerical results and discussion 

In the following we shall use reduced quantities: p* = pa^ as the reduced density, 
m* = ma/ y/ksTa^ as the dimensionless dipole moment of species a, H* = H^/ a'^ / (ksT) 
as the reduced magnetic field strength, and M* = / (ksT) as the reduced 

magnetization. The calculation of the zero-field susceptibility and the magnetization 
M{H) of the multicomponent DHS fluid mixtures (with identical particle diameters but 
different dipole moments) can be summarized by the sequence of the following steps: 

1) calculation of the Langevin susceptibility according to Eq. Qj, 

2) solving Eq. §i) for f, 

3) calculation of the zero-field susceptibility according to Eq. (|TT]l . 

4) calculation of the magnetization for a given value of H according to Eq. (129|) using 
the consecutive approximation method with M = as initial value. The convergence of 
this consecutive approximation is very good, obtaining the limiting results within 5-8 
cycles. 

Figure 1 shows the dependence of the zero-field susceptibility x the Langevin 
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Figure 3. Same system as in Fig. 2. For xi = 0.25 and xi = 0.75 DFT (MSA) and 
the corresponding MC data are compared witlr tlie Langevin magnetization curves. 
The error bars of the MC data are given by the symbol sizes. 



susceptibility Xl ~ (^^^ -^l- ©) dipolar hard sphere mixtures as obtained 
from Eq. (fTTIl and from the numerical solution of Eq. ([8]). This result is compared 
with MC simulation data for a binary mixture [{m\,m2) = (0.5,1)] for three total 
number densities p* and five concentrations Xi. In these cases DFT (MSA) provides 
a good approximation for the initial susceptibility of this binary system within the 
range < 47rx^ < 2.5. Within this range the two-component system with various 
concentrations and densities can be described by the same master curve which is the 
same also for different systems {ml, rn^j). This is the main statement of the MSA 

theory by Adelman and Deutch |23] . 

Figure 2 displays the magnetization curves of the two-component DHS fluid mixture 
with {ml, 7712, p*) — (O-S, 1,0.6) for five values of the concentration. For high values of 
H* we find excellent, quantitative agreement for all concentrations between the DFT 
(MSA) results and the MC data. For small values of H*, especially for H* = 0.5, i.e., 
in the linear response regime, the agreement between the simulation data and the DFT 
results is also very good, which matches with the good agreement found for the zero- 
field susceptibility (see Fig. 1). Close to the elbow of the magnetization curves the 
level of quantitative agreement reduces significantly for smaller concentrations Xi of the 
magnetically weaker component, while it remains good for large concentrations xi. We 
note that also for two-dimensional systems this range is the most sensitive one concerning 
the agreement between theoretical results and simulation data |26] . For the same system 
Fig. 3 displays a comparison between the DFT results together with the MC data 
and the corresponding Langevin theory. This shows that the interparticle interaction 
enhances the magnetization relative to the corresponding values of the Langevin theory. 
For a three-component DHS fluid mixture [{m^l, mg, m^) = (0.5, 0.75, 1)] Fig. 4 displays 
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Figure 4. Dependence of the magnetization on the concentration 2:3 of the species 



with '1 



1 for a three-component DHS fluid mixture for a fixed field strength H* 



and concentrations xi = 0.25, xi = 0.5, and xi — 0.75; xi + X2 + X3 = 1. The reduced 
dipole moments of the species 1 and 2 are m| = 0.5 and — 0.75. The overall 
number density is p* = 0.6. The error bars of the MC data are given by the symbol 
sizes. 



the dependence of the magnetization on the concentration X3 for a fixed field strength 
H* = 2 and for three values of the concentration xi, Xi + X2 + = 1. Since the 
value H* = 2 falls into the aforementioned elbow regime of the magnetization curves, 
for the small concentration xi = 0.25 of the less polar {ml = 0.5) component the DFT 
(MSA) results underestimate the simulation data. At higher concentrations of Xi the 
agreement between DFT and the simulation data is much better. This is expected to 
occur, because for large Xi the fluid is dominated by less polar particles. 

6. Summary 

We have obtained the following main results: 

1) Based on a second-order Taylor series expansion of the anisotropic free energy 
functional of equally sized dipolar hard spheres with different dipole moments and by 
using the mean spherical approximation (MSA) we have derived an analytical expression 
(Eq. (B9]) ) for the magnetization of multicomponent ferrofluidic mixtures in external 
fields. This implicit equation extends the applicability of MSA to the presence of 
external magnetic fields of arbitrary strengths. We find quantitative agreement between 
the results from this DFT (MSA) and our Monte Carlo simulation data for Langevin 
susceptibilities 4:TcXl ^2.5 (Figs. 2, 3, and 4). 

2) As confirmed also by MC simulation data the zero-field susceptibility of 
multicomponent ferrofluids can be expressed by a single master curve in terms of the 
Langevin susceptibility (Fig. 1). Beyond the linear response regime the magnetization 
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curves of multicomponent ferrofluids cannot be reduced to a single master curve (Eq. 



3) By applying the MSA theory for the magnetic susceptibility to polydisperse systems 
we have extended the multicomponent magnetization equation to polydisperse systems. 
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